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We present a dual-resolution molecular dynamics (MD) simulation of liquid water employing a 
recently introduced Adaptive Resolution Scheme (AdResS). The spatially adaptive molecular resolu- 
tion procedure allows for changing from a coarse-grained to an all-atom representation and vice-versa 
on-the-fly. In order to find the most appropriate coarse-grained water model to be employed with 
AdResS we first study the accuracy of different coarse-grained water models in reproducing the 
structural properties of the all-atom system. Typically coarse-grained molecular models have a 
higher diffusion constant than the corresponding all-atom models due to the reduction in degrees of 
freedom (DOFs) upon coarse-graining that eliminates the fluctuating forces associated with those 
integrated-out molecular DOFs. Here, we introduce the methodology to obtain the same diffusional 
dynamics across different resolutions. We show that this approach leads to the correct description of 
essential thermodynamic, structural and dynamical properties of liquid water at ambient conditions. 



I. INTRODUCTION 

Water is a typical example where the interplay between 
different length scales determines the relevant properties 
of the system. Including explicit water in large biomolec- 
ular simulations is crucial^ but normally not feasible 
due to the large number of water molecules needed to 
describe biomolecular function. The computational re- 
quirements of all-atom simulations in explicit water usu- 
ally do not allow to investigate biologically relevant time- 
scales^. Most computational approaches renormalize the 
role of water into the definition of effective inter-residue 
interactions^ or as a continuous field, in order to focus the 
computational efforts to simulate the biomoleculc. How- 
ever, the discrete nature of water molecules can play a 
role in, e.g., controlling protein functionality or the case 
of interactions between hydrocolloids and water, where 
hydrogen bonding is important, etc. To overcome this 
problem, a multiscale modeling of water can be very ad- 
vantageous to speed up the simulation without losing 
physics-chemical accuracy. 

Multiscale modeling is emerging as a viable way 
to bridge the various length and time scales involved 
in complex molecular systems^^&^ i 10 ' 11 ' 12 ' 13 ' 14 . In 
general, to build a multiscale model two main issues 
need to be addressed. The first issue consists in mapping 
the system into a robust reduced model while preserving 
the here relevant physico-chemical properties, i.e., radial 
distribution functions, pressure, and temperature, of the 
reference all-atom system. The second issue involves the 
definition of a robust and physically accurate procedure 
to smoothly join the different resolutions. We have 
recently proposed the Adaptive Resolution Scheme 
(AdResS) 15 ' 16 to address both issues for a system of 
water molecules^. Typically, the reduction in the 
number of the system degrees of freedom (DOFs) upon 
coarse-graining introduces a time scale difference be- 
tween the coarse-grained and explicit region. Although 



in certain instances the difference in time scale may be 
advantageous for reaching longer simulation times, in 
other cases having the correct diffusional dynamics is 
crucial. This is the case, for instance, if a dynamical 
property is focus of an investigation as in the transloca- 
tion of biopolymers through membrane nanopores 18 . 

In this paper, in order to justify our choice of the 
recently introduced single-site water modelii as the 
optimal coarse-grained water model to be used with 
AdResS (when using the standard three-site TIP3P wa- 
ter models as the all-atom water model), we present 
firstly a detailed analysis and comparison of newly de- 
veloped coarse-grained water models having a different 
number of DOFs (from 3 to 9) with the single-site wa- 
ter modelii as well as the TIP3P water models. In this 
way we can study the contribution of separate DOFs of 
the all-atom water model to the definition of its struc- 
tural properties. As it turns out the single-site water 
model can reproduce the structural and thermodynami- 
cal properties of the all-atom system that are relevant for 
the multiresolution simulation equally well as the more 
sophisticated models. Therefore, we have chosen it as the 
coarse-grained model in our adaptive resolution simula- 
tions. Next, we introduce the methodology to obtain the 
same diffusional dynamics across the different resolutions 
in the hybrid atomistic/coarse-grained (ex-cg) model sys- 
tem composed of explicit and coarse-grained molecules as 
presented in Figure [TJ 



II. MAPPING OF A COARSE-GRAINED 
MODEL TO THE ALL-ATOM 
REPRESENTATION 

In order to map a coarse-grained to an all-atom model, 
we have to construct an effective potential between 
coarse-grained molecules in a way that the atomistic 
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FIG. 1: (Color online) Top figure: The coarse-grained molecule is rep- 
resented at the left, and the all-atom water molecule is represented at 
the right. The middle hybrid molecule interpolates between the two. 
Bottom figure: The dashed line represents the weighting function w(x) 
G [0, 1] defined in Ref^— and discussed in the text. The values w=l 
and w— correspond to the atomistic and coarse-grained regions of the 
hybrid atomistic/mesoscopic system. The full line represents the inter- 
face correction weighting function s(w) denned in Figurc[6] The value 
s=l corresponds to the atomistic and coarse-grained regions, while s— 
when w— 1/2. The explicit and interface regions arc divided in three 
regions ( A, B , C) of width 4 A for computing the cosine distribution in 
Figure [8] 



structural properties are reproduced. We numerically 
build such effective Hamiltonian following an "inverse 
statistical-mechanics" approac h 20 ! 21 1 22 . In order not to 
bias a priori the potential function with the choice of a 
particular functional form, we introduce a grid approxi- 
mation, expressing a general pairwise Hamiltonian in the 
form: 



H — ^ U a S a 



(i) 



where S a is the number of coarse-grained particles in 
the interval € (r a , r a + dr) and U a is the discretized 
value of the potential in the interval (r Q ,r Q + dr). The 
averages (S a ), which can be interpreted as a discretized 
radial distribution function (rdf), are functions of the set 
of parameters {U a }. Values of the potential parameters 
{U a } that reproduce the rdf of the all-atom model can 
be obtained iteratively as follows: 

• The center-of-mass rdf g(r) is computed from all- 
atom simulations and used as a target function. 
The potential of mean forces is assumed as an initial 
approximation to the effective potential function: 



U° a = -kT\ng(r a ) 
where U a is the potential at a distance r D 



(2) 



By comparing the atomistic center-of-mass rdf with 
the rdf obtained for the coarse-grained model at 
the n-th iteration (with corresponding potential pa- 
rameters f/™) the (n+ l)-th correction to the effec- 
tive potential can be found. Namely, the correction 



to the potential parameters {/" is given by the so- 
lution of the system of linear equations^: 



A{S a )=J2 



where: 

d(S a ) 



9Kv 



= -/? (S a S 7 )-(S a )(S 7 ) 



(3) 



(4) 



Ensemble averages (S a S*y) and (S a ),(S-y) are com- 
puted from MD simulations (with potential param- 
eters [/"), and the correction A[7™ to the potential 
£7™ is obtained by solving the system of linear equa- 
tions 



jjn+l V n , Au r, 



(5) 



• The corrected potential parameters U™ +1 are then 
used to perform molecular dynamics simulation 
with the coarse-grained model and calculate new 
ensemble averages for the quantities (S a S 7 ) and 

(S a ) ,(S-y) . 

The points above are repeated to determine a new set 
of corrections At/™ +1 . The procedure is repeated until 
convergence is reached within the statistical simulation 
error. 

Additionally, to match the pressure of the coarse-grained 
to the all-atom model, after each iteration, a weak con- 
stant force is added to the effective force in such a way 
that the total effective force and potential energy are zero 
at the fixed cutoff distance r, i 16 i 23 . 



AU(r) = U [1 



(0) 



We use a value of r c = 7 A as the rdf 1 for r > r c . 
Depending on the pressure in the current iteration be- 
ing above or below the target value corresponding to the 
pressure of the reference all-atom system, U Q can assume 
a positive or negative value. 

For the atomistic representation of water we selected the 
rigid TIP3P models that gives a reasonable description 
of the behavior of liquid water around standard (phys- 
iological) temperature and pressure, and is not too ex- 
pensive computationally. The potential function of rigid 
TIP3P involves a rigid water with three interaction sites 
that coincide with the atomic positions. The model 
uses atom-centered point charges to represent the elec- 
trostatic, i.e., positive charges on the hydrogens and a 
negative charge on oxygen (eo — — 2e#). The van der 
Waals interaction between two water molecules is mod- 
eled by the Lennard-Jones potential with just a single in- 
teraction site per molecule centred on the oxygen atom; 
the van der Waals interactions involving the hydrogens 
are omitted. This gives the following potential: 
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The values of the parameters e, A, and C are assigned to 
reproduce reasonable structural and energetics properties 
of liquid water—. 



III. COARSE-GRAINED MODELS FOR WATER 

Many simplified coarse-grained models have been de- 
veloped to reproduce qualitatively the structural prop- 
erties of water—. Existing coarse-grained models can 
of course only reproduce thermodynamical properties of 
all-atom water to a certain exten d 25 ' 26 i 27 ' 28 i 29 . By us- 
ing the procedure described above, we have designed wa- 
ter models at different level of coarse-graining (three-site, 
two-site, and one-site) schematically depicted in Figure 
[2j The gradually increasing level of coarse-graining al- 
lows us then to study the contribution of different DOFs 
to the structural and thermodynamical properties of the 
reference all-atom system. In this way we can determine 
the most appropriate coarse-grained model to be used in 
our adaptive resolution simulations. 



penalty function defined as 




FIG. 2: (Color online) Cartoon of different water models employed in 
this study. From top to bottom: The rigid TIP3P water model, and 
three-site, two-site, one-site coarse-grained water models, respectively. 



Three-site interaction model 

A first level of coarse-graining can be introduced by 
preserving the atomic positions of each atom in the 
TIP3P water model and replacing the explicit electro- 
static interactions (the first term in Eq. ([7])) with ef- 
fective short-ranged interactions. This model preserves 
explicitly the H-bond directionality but it is computa- 
tionally less expensive than TIP3P because it removes 
the long-ranged Coulombic term. The calculated all- 
atom site-site rdfs are used as input to build the 00, 
OH and HH effective potentials (see Eqns. [4]). In order 
to quantify the agreement of the rdfs (corresponding to 
the all-atom and coarse-grained model) we introduce the 



fp 



9 C9 (r/aoo)-9 ex (r/aoo) 



exp(-r/a o)d(r/a o), 
(8) 



where g ex and g cg are the reference site-site rdf of the 
atomistic and coarse-grained system, respectively, and 
(Too is the Lennard-Jones constant of the TIP3P water 
model^. The extremely low values of the f p obtained 
(f p of the 00 rdf is 7.9 * 1(T 5 , for OH is 8.2 * 1CT 5 and 
HH is 1.1 * 1CP 4 ) for the optimized effective potential for 
the three-site model indicates a perfect matching of the 
rdfs. The optimized effective potentials are as expected 
very similar to the effective potentials of other previously 
introduced 3-site coarse-grained water model o 28 i 30 . In 
addition, to check the angular properties (that are not 
completely defined by the rdf) we computed the distri- 
bution of the angle 6 formed between the center-of-mass 
of three nearest neighbor molecules (Figure [3J top) , and 
the distribution of the orientational order parameter q 
(see Figure UJ), as defined by Errington et al~: 



o 3 4 
j=l k=j+l 



COSIpjk + - 



(9) 



where ipjf. is the angle formed by the lines joining the 
oxygen atom of a given molecule and those of its nearest 
neighbors j and k. The parameter q measures the extent 
to which a molecule and its four nearest neighbors adopt 
a tetrahedral arrangement. The excellent agreement be- 
tween the rdfs and the angular distribution suggests that 
the explicit H-bond directionality is responsible for the 
local structure. On the contrary, the small deviation of 
the orientational order parameter distribution suggests 
that the electrostatic interactions are at least partly re- 
sponsible for the overall tetrahedral arrangement of the 
system. 



Two-site interaction model 

A two-site simplified model of water changes the sym- 
metry of the molecule by using two interacting sites in- 
stead of three, one corresponding to the oxygen atom and 
the second one providing the directionality of the dipole 
moment. Since the dipole moment of TIP3P is [i = 2.35D 
and the charge separation is q=0.82ei£, the effective sep- 
aration of the negative and positive charge centers is 
d = & re 0.5966 A. Thus, we define the two-site model by 
representing each water molecule by two spheres, which 
are joined by a rigid rod of length 0.5966 A. Two-site 
models have been previously proposed for water—, how- 
ever existing models cannot exactly reproduce the radial 
structure of the reference all-atom water model. The two- 
site water model with parameters optimized by the proce- 
dure described above can reproduce remarkably well the 
rdfs for the different pairs of interactions: oxygen-oxygen 
with f p = 4.5 * 10~ 5 , oxygen-dipole with /„ = 1.4 * 10~ 4 
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and dipole-dipole with f p = 3.2 * 10~ 5 . Figure [3] shows 
the angular distribution and the orientational order pa- 
rameter distribution obtained for the two-site model, in 
comparison with the results of the all-atom and three- 
site models. The excellent agreement between the rdfs 
suggests that the inclusion of a dipole center is sufficient 
for the correct radial structure. The deviation on the 
angular distribution from the three-site model indicates 
that the explicit H-bond directionality has a contribution 
to the angular structure (in addition to the electrostatic 
contribution, as found in the three-site model). 



One-site interaction model 

In the one-site coarse-grained model, the water is rep- 
resented by one spherically symmetrical site having a 
mass mo + 2m h. As demonstrated in Refji 7 -, the one- 
site model can reproduce the center-of-mass rdf (f p — 
7.6* 10~ 5 ) and thermodynamic properties of the reference 
all-atom water model with remarkable accuracy. The op- 
timized effective potential (see inset of Figure [5]) has a 
first primary minimum at about 2.8 A corresponding 
to the first peak in the center-of-mass rdf. A second, 
slightly weaker and significantly broader minimum at 4.5 
A corresponds to the second hydration shell. The com- 
bined effect of the two leads to a local packing close to 
that of the all-atom TIP3P water. Our effective coarse- 
grained potential is quite different from the previously 
suggested potential o 25 i 26 i 27 i 31 : while in previous one-site 
models the deepest minimum corresponds to the second 
hydration shell, the absolute minimum in our model is 
found in the first shell. 

The good structural agreement between the explicit 
and coarse-grained models indicates that although our 
coarse-grained model of water is spherically symmetric 
and therefore does not have any explicit directionality, it 
approximately captures the correct local structure. Fur- 
thermore, the angular and orientational order parameter 
distribution coincides with what obtained for the two-site 
model. Since the one-site model can qualitatively predict 
the correct angular distribution and is computationally 
the least expensive of all the coarse-grained models, it 
was our choice for the multiscale approach (presented 
mil). 

As a consequence of the reduced number of DOFs, 
there is a time scale difference in the dynamics of the 
coarse-grained system with respect to the atomistic one. 
The diffusion coefficient increases as the level of coarse- 
graining is increased, with the diffusion coefficient of the 
center of mass for the TIP3P model D « 3.7 * 1CT 5 ^, 

S ' 

_ r 2 

for the three-site model D w 4.68* 10~ o£21 ^-, for the two- 

s " 

r 2 

site model D ^ 7.54* 10~°^— and for the one-site model 

s 

D « 8.1* 1CT 5 ^. 

s 

Using the one-site coarse-grained model one can sig- 
nificantly speed up the simulation, with a total gain in 
computational time of a factor ~ 17 — 20 when compared 
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FIG. 3: Top figure: The center-of-mass angular distribution between 
three nearest neighbors for all-atom, 1-site coarse-grained, 2-site coarse- 
grained and 3-sitc coarse-grained models for water. Bottom figure: The 
analogous distributions of the orientational order parameter q. The 
center-of-mass rdf of the different water models coincides exactly with 
the all-atom rdf (not shown) . 

to atomistic simulations 17 . This is due to the reduction 
of the number of interactions, which are also softer than 
in the atomistic case, and due to an intrinsic time scale 
difference in the diffusive dynamics of the coarse-grained 
system, that is faster of about a factor 2 than the all- 
atom simulations counterpart (see the above values for 
the diffusion constant). 



IV. MATCHING THE DIFFUSIVE DYNAMICS 
OF THE COARSE-GRAINED TO THE 
ALL-ATOM MODEL 

As discussed above, the dynamics of the one-site 
coarse-grained model is faster than that obtained from 
all-atom simulations. The speed-up occurs because of 
the reduction in DOFs upon coarse-graining, which elimi- 
nates the fluctuating forces associated with those missing 
molecular DOFs, leading to the much smoother overall 
energy landscap e) 33 ! 34 . In our MD simulations we use a 
Langevin thermostat, and the equations of motion are in 
the form: 



CLV ' 

m~ = F t - rrnVvi + R l (t) (10) 

where Fi — — fj 1 is the deterministic force and 
Ri is a stochastic variable with (Ri(t)Rj(t + r)) = 
2TmikTS(T)Sij. The coefficient T determines the 
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strength of the coupling to the bath [not to be confused 
with the friction coefficient £ of the system (for TIP3P 
£ = 288. 6ps -1 )]. When jt is large compared to the typ- 
ical time scales in the system, the stochastic dynamics 
reproduces dynamics obtained by molecular dynamics 
(MD)2£, while if ^ is small, then the stochastic dynamics 
deviates significantly from MD. Simulations performed 
using the the one-site model with same coupling to the 
bath (r = bps- 1 < that is usually used for TIP3P 36 
produce an increase of a factor of about 2 in the diffu- 
sion time scale of the coarse-grained system. An acceler- 
ated dynamics can be advantageous in some cases but as 
mentioned previously it can be a problem if the dynam- 
ical properties themselves are under investigation. The 
coarse-grained dynamics can be slowed down by increas- 
ing the effective friction in the coarse-grained system. 
That is, by changing T, the one-site model can mimic 
the diffusive dynamic of the all-atom system. Figure [4] 
shows that the diffusion coefficient D monotonically de- 
creases as the coefficient T increases. The coarse-grained 
model yields the same diffusion coefficient as the all-atom 
system (for which T ex — bps -1 ) when T cg — lbps -1 . The 
additional frictional noise does not affect the short time 
dynamics of the system since T ex < T cg <C It is 

worth mentioning that the structural/thermodynamics 
of the coarse-grained system are not effected by changes 
in r in this range. Figure 2] shows the perfect agreement 
between the center-of-mass rdf for the coarse-grained and 
all-atom system when different values of T are used for 
the two models. However, we have to bear in mind that 
by employing the Langevin thermostat, Eq. ([1(7)) . we 
screen the hydrodynamics. Furthermore, for too large T 
the statics is also altered because one arrives eventually 
at the Brownian dynamics. In order to correctly repro- 
duce the hydrodynamics one has to resort to the Dissipa- 
tive Particle Dynamics (DPD) thermostat as was done 
for instance in Ref. — for the case of a macromolecule in 
the hybrid atomistic/mesoscale solvent. 



V. ADAPTIVE RESOLUTION SCHEME 
(ADRESS) 

As we have proposed ir*i£, we design an adaptive multi- 
scale system where half of the simulation box is occupied 
by atomistic water molecules and the other half by the 
corresponding one-site coarse-grained molecules, respec- 
tively, as schematically presented in Figure [TJ In order to 
smoothly couple the regimes of high and low level of de- 
tail, we apply the AdResS schem o 15 ' 16 ' 17 , that allows the 
molecules to move freely between regimes without feeling 
any barrier in between. The interface region contains hy- 
brid molecules that are composed of an all-atom molecule 
with an additional massless center-of-mass particle serv- 
ing as an interaction site. The transition is governed by a 
weighting function w that interpolates the molecular in- 
teraction forces between the two regimes, and assigns the 
identity of the molecules. In the present work we used the 
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FIG. 4: The diffusion coefficient as a function of the friction constant. 
The triangle indicates the value obtained for the all-atom diffusion co- 
efficient when a value T — 5ps — 1 is used in the simulation. The black 
dots indicate the diffusion coefficient obtained for the coarse-grained 
system for different values of T. The two circles indicate the diffusion 
coefficients obtained for the coarse-grained system when T — 5ps~ 1 and 
r — 15ps -1 , respectively, are used in the simulation. The inset com- 
pares the rdfs obtained for the all-atom system and the coarse-grained 
system for different values of the friction coefficients. 



weighting function defined in^ and described in Figure 
[TJ The function w is defined in such a way that w = 1 
corresponds to the atomistic region, and w — to the 
coarse-grained region, whereas the values < w < 1 cor- 
respond to the interface layer. This interpolation leads 
to intermolecular forces acting on the center-of-mass of 
the molecules a and as: 

F a/3 = w(X a )w(X )F^ m + [1 - w(X a )w(Xp)]F%S. 

(11) 

F Q/ 3 is the total intermolecular force acting between cen- 
ter of mass of molecules a and (5. F^g m is the sum of all 
pair atomic interactions between explicit water atoms of 
molecule a and explicit water atoms of molecule /?, and 



a(3 



is the effective pair force between the center-of-mass 
of two water molecules. X a and Xp are the center-of- 
mass coordinates of molecules a and (3. Note that the 
AdResS as given by Eq. ([TTj) satisfies Newton's Third 
Law. This is crucial for the diffusion of molecules across 
the resolution boundaries 1 ^. Since the total force as de- 
fined by Eq. ([TTj) depends on the absolute positions of 
the particles and not only on their relative distances it is 
in general not conservative, i.e., the work done by it de- 
pends on the path taken in the transition regime. Hence 
the corresponding potential does not exist and the total 
potential energy of the system can not be define d 39 ' 40 . 
Still, the intermolecular forces between molecules outside 
the transition regime are conservative with well defined 
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potentials, i.e., the all-atom or effective coarse-grained 
potentials^ 5 .. 

Each time a molecule crosses the boundary between the 
different regimes it gains or looses (depending on whether 
it leaves or enters the coarse-grained region) its equili- 
brated rotational DOFs while retaining its linear momen- 
tum. By this choice of interactions in the interface re- 
gion, the hybrid molecule interacts with molecules in the 
coarse-grained region on a coarse-grained level. On the 
other hand, the interactions of the hybrid molecules with 
the molecules in the explicit region are a combination of 
the explicit-explicit and cg-cg interactions to smoothly 
and efficiently equilibrate additional DOFs upon mov- 
ing in the explicit regim o 15 ' 16 . This change of resolution, 
which can be thought in terms of similarity to a phase 
transition^, requires to supply or remove latent heat and 
thus must be employed together with i.e. a Langevin 
thermostatic. The thermostat is coupled locally to the 
particle motion and provides a mean to deliver or absorb 
the latent heat. 

As in Reffii the reaction field (RF) method is used, 
in which all molecules outside a spherical cavity of a 
molecular based cutoff radius R c = 9 A are treated as 
a dielectric continuum with a dielectric constant €rf = 

^,43,44,45,46,47 The Coulomb force acting on a part ial 

charge ei a , belonging to the explicit or hybrid molecule 
a, at the center of the cutoff sphere, due to a partial 
charge ej , belonging to the explicit or hybrid molecule 
(3, within the cavity is: 



•c ic 



(J, 
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x3p- 



To suppress the unphysical pressure fluctuations 
emerging as artifacts of the scheme given in Eq. [TT] we 
employ an interface pressure correction— within the tran- 
sition zone. The latter requires a re-parametrization of 
the effective potential in the system composed of exclu- 
sively hybrid molecules (with a constant of w=l/2). We 
redefine the center-of-mass intermolecular forces as 

f$ = s K^)^)]p^+(i-«[t t ftWx <) ])jq (t) 

where the function s G [0,1] is defined as 

J COs(7Ty / x) 2 



if < x < 0.25 
if 0.25 < x< 1.0 



(14) 



so that s[0]=l, s[l]=l, and s[l/4]=0, as shown in Figure 
The interface correction force F^a. is the effective pair 
force between molecules a and defined by the effective 
pair potential, which is obtained by mapping the hybrid 
model system composed of exclusively hybrid molecules 
with w — 1/2 to the explicit model system. The corrected 
effective potential J7f c m shown in Figure [5] is obtained by 
mapping the ex — cg(w — 1/2) j c system containing only 
hybrid molecules. The minima of the effective poten- 
tial Uf™ become deeper, and a higher barrier separates 



the first and second hydration shell when compared to 
the center-of-mass effective potential U cm . The center- 
of-mass rdf of the hydrid system exactly reproduces the 
structure of the reference system as shown in Figure O 
The mixing function in Eq. [14] can correctly reproduce 
the thermodynamic properties of the (ex — eg) system 
with the same mean temperature (0.1 % of difference) 
and pressure (0.5 % of difference). Detailed comparisons 
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FIG. 5: The center-of-mass rdfs for explicit, cx-cg(w— 1/2), and 
coarse-grained systems. The inset shows the corrected effective pair 
potential for hybrid molecules [dashed line] and the reference ef- 

fective potential for the coarse-grained molecules [full line]. 



between the bulk explicit simulations and the explicit 
regime in our hybrid setup prove that this approach docs 
not alter the structural properties of the water model 
studied. In particular, Figure [7] shows that the struc- 
tural properties of the explicit regime in the multiscale 
system are exactly the same as in bulk explicit simula- 
tions. No change is detected in the orientational prefer- 
ences of water molecules near the interface region. It is 
worth mentioning that similar results are obtained in the 
interface region (results not shown). 

Figure [5] shows the probability density function of the 
orientational parameters cos 8oh (O-H vector and the 
normal of the interface surface) and cos Jhh (H-H vector 
and the normal of the interface surface) in three consec- 
utive layers of width 4 A next to the interface region 
along the X-axis (see Figure [1] for the definition of the 
regions). The uniform cosine distribution of the water 
molecules in layers inside and next to the interface indi- 
cates that the orientational DOFs are fully equilibrated 
in the hybrid region, as the molecules do not have to re- 
orient upon crossing the interface region and there is no 
change of behaviour at the interface^. 
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FIG. 6: The interface correction weighting function s(w). The value 
s— 1 corresponds to the atomistic and coarse-grained regions and s— 
when w— 1/2. 




FIG. 7: The ccntcr-of-mass (fp=7. 65 * 10 -5 ), OH (fp=6 * 10~ 5 ) and 
HH (fp— 6.6 * 10~ 5 ) rdfs for the explicit region in the hybrid system 
[dots], and bulk [line] systems. 



VI. POSITION DEPENDENT THERMOSTAT 

As mentioned above, the structural/thermodynamics 
of the coarse-grained system are not changed by increas- 
ing the background friction in the Langevin thermostat. 
To obtain the same diffusional dynamics across different 



FIG. 8: (Color online) Cosine distribution of the angle formed by the 
H-H vector (top panel) and O-H vector (bottom panel) of the water 
molecules with the interface normal vector pointing toward the explicit 
region in three consecutives regions of width 4 A each (see Figure 
for the definition of the regions). 



resolutions, the coefficient T in the Langevin thermostat 
is changed on-the-fly depending on the number of DOFs 
of the molecules. As shown in Figure^ when a constant 
coefficient T is used, two regimes are observed in systems 
composed of only hybrid molecules: 

• For w e [0 — 0.6] the value of the diffusion coefficient 
D is essentially constant; 

• For w e ]0.6 — 1.0] the value of the diffusion coeffi- 
cient D drops steeply with w. 

When the molecular identity w is greater than 0.6 the 
hybrid molecules start to equilibrate their orientational 
structure and their diffusive dynamics is slowed down, 
as shown in Figures [9] and [TDJ In order to obtain the 
same diffusive dynamics for different levels of resolution 
we propose the following functional form for the friction 
coefficient in the Langevin thermostat: 



F C9 if.<0.6 
v ' \aw + P if 0.6 < w< 1.0 



(15) 



This choice provides a simple interpolation between the 
two limit values of r(0.6) = T(0) = T cg = 15ps _1 and 
r(l) = Tan-atom — 5ps . The parameters a and (3 are 
— 25ps~ 1 and 30ps _1 , respectively. As shown in Figured 
systems containing only hybrid molecules with different 
particle identities exhibit the same diffusional dynamics 
when T(w) has the functional form proposed in Eq. 1151 
In our simulation we employ a local Langevin thermo- 
stat with a position dependent frictio n 48 i 49 i 50 i 51 to match 
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FIG. 9: Top figure: The dashed curve indicates the dependency of 
the friction coefficient as a function of the particle identity w when a 
position dependent thermostat is used. The full line shows the constant 
value of the friction coefficient when a regular thermostat is used. Bot- 
tom figure: The dots indicate the diffusion of the molecules when the 
regular thermostat is used. The triangles indicate the diffusion of the 
molecules when the position dependent thermostat is used. 



the diffusion constants of the coarse-grained and all-atom 
regimes. The Langevin equation with a position depen- 
dent coefficient T(x) can be written as^: 

rriidvi/dt = Fi — m l T(x)vi + Ri(x,t) (16) 

where Ri(x, t) is: 

(Ri(x,t))=0, (17) 

(Ri(x,ti)Rj(x,h)) = 2T{x)m l kT5{t l - t 2 )5 lJ (18) 

The functional form of T(x) depends on the weighting 
function w(x) and it is shown in Figure [9] (see Eq [T5|) . 
There is no difference in the structural and thermody- 
namic properties of the system when a position depen- 
dent T(x) is used instead of a constant coefficient T. The 
density is homogeneous in both the coarse-grained and 
explicit regions with (very) small oscillations in the tran- 
sition regime, cf. Figure [TT1 . 

To prove the free exchange of molecules between the dif- 
ferent regimes we have computed the time evolution of a 
diffusion profile for molecules that were initially localized 
within two slabs of width 9.5 A neighboring the in- 
terface layer. The molecules initially localized within the 
two slabs spread out symmetrically with time when a po- 
sition dependent friction is used in the local thermostat. 
On the contrary, the molecules spread out asymmetri- 
cally with time as shown in Figure [12] when a fixed coef- 
ficient r is used. This asymmetry arises from the above- 
mentioned difference in diffusion coefficient D between 
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FIG. 10: (Color online) The H-H.O-H and O-O rdfs for different par- 
ticles identities. The water molecules start to correct their structural 
properties when w > 0.6. 



the all- atom and coarse-grained regions (see Ref>^). Fig- 
ure [13] shows the behaviour of the diffusion coefficient D 
as a function of the x position using particles that are 
within a slab of « 9.5 A of the mixed system for the two 
different type of thermostats. As Figure Q1J] illustrates, 
the molecules are slowed down when they cross the in- 
terface from the coarse-grained region to the explicit re- 
gion if a constant friction is used for the thermostat. The 
change in the diffusion coefficient D is localized in the in- 
terface region while D is approximately constant inside 
each region, coarse-grained [D « (8.04 ±0.54) x 10 -92 ^-] 

and explicit [D « (4.2 ± 0.28) x 10~ 9 ^]. When a po- 
sition dependent friction is used for the local thermostat 
the diffusive dynamics of the molecules is the same as for 
the all-atom system across all regions, as shown in Figure 



VII. CONCLUSIONS 

In this paper we extended the multiscale model for 
water we have recently proposed^ to obtain the same 
diffusional dynamics across the different resolutions. We 
also studied the accuracy of different coarse-grained wa- 
ter models in reproducing the structural properties of 
the all-atom system. The results of this study show that 
for our purposes the single-site model performs equally 
well as the more sophisticated (three-site and two-site) 
coarse-grained models and is hence the most appropri- 
ate for the presented adaptive resolution simulations. 
We envisage that such adaptive resolution simulations 
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FIG. 11: Top figure: Normalized density profile in the x-dircction of 
the mixed system when the position dependent thermostat is used. Bot- 
tom figure: Normalized density profile in the x-direction of the mixed 
system when a constant friction constant is used. 



of water will play an important role in the modeling of 
biomolecular systems where the coupling of different time 
and length scales is crucial to understand their physico- 
chemical properties. 
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VIII. APPENDIX 

A. Multiscale simulation protocol 

All the results presented in the paper were obtained by 
performing nVT simulations using the ESPResSc 52 (for 
the one-site and multi-scale model) and AMBER 8— (for 
the 2-site and 3-site models) simulation packages with a 
Langevin thermostat, with a friction constant r = 5 ps -1 
when a regular thermostat is used and a time step of 
0.002 ps. All the models considered a system of of 1464 
water molecules at T re f — 300 K and p = 0.96 g/cm 3 . 
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FIG. 12: (a) Time evolution of a diffusion profile when the position 
dependent thermostat is used for molecules that arc initially (at time 
t — ps) located inside two neighboring slabs at opposite sides of the 
mid interface layer. The diffusion profile is averaged over ~ 400 differ- 
ent time origins. Vertical lines denote the boundaries of the interface 
layer. Bottom figure: The diffusion profile for the all-atom side is av- 
eraged over S3 400 different time origins. Middle figure: The same as 
before but for molecules that are initially localized inside the slab on 
the coarse-grained side of the interface region. Top figure: Difference 
between the corresponding diffusion-profile of the coarse-grained and 
the all-atom regions, (b) Time evolution of a diffusion profile with the 
regular thermostat using the same criteria as given in (a). 



The density was obtained from an all-atom NPT simula- 
tion with P re f = 1 atm using a Reaction-field method for 
the electrostatics and a cut-off method for the pressure. 
The results presented for the 2-site and 3-site models 
were obtained for a density p = 1.007 g/cm 3 . This value 
of the density was obtained from an all-atom NPT simu- 
lation with P re f = 1 atm using an Ewald-method for the 
electrostatics and a long-ranged correction for the pres- 
sure after the cut-off. Periodic boundary conditions and 
minimum image convention were applied in all directions. 
The bonds and angle of the water molecules were con- 
strained by using the RATTLE procedure. After warm- 
up and equilibration, several trajectories of 1.5 ns each 
were collected for each different model. All simulations 
were performed with a force capping to prevent possible 
force singularities that could emerge because of overlaps 
with neighboring molecules when a given molecule enters 
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FIG. 13; Diffusion coefficient as a function of the x position using 
particles that are within a slab of ~ 9.5 A of the hybrid system for 
different type of thermostats. Each point in the plot is centered in each 
slab. The average is over ~ 400 different time origins. The horizontal 
line is the diffusion coefficient of a system containing only all-atom 
water molecules. The dots indicate the diffusion of the molecules across 
the hybrid system when the regular thermostat is used. The triangles 
indicate the diffusion of the molecules across the hybrid system when 
the position dependent thermostat is used. 



the interface layer from the coarse-grained side. The size 
of the system is 94.5 A in the x direction and 22 A in 
both the y and z directions. An interface layer of 18.9 A 
between the coarse-grained and all-atom models is set 
along the x direction. To treat all the molecules of the 
system equally regardless of their level of detail, we define 
the pressure of the system according to their lower level 
of detail, that is by the molecular pressur o 15 ' 16 > 54 ' 55 . The 
temperature was evaluated using the fractional analog of 
the equipartition theorem: 



(K a ) = 



aksT 



(19) 



where (K a ) is the average kinetic energy per fractional 
quadratic DOF with the weight w(r) = a^-. Via Eq. ([!§)) 
the temperature is also rigorously defined in the transi- 
tion regime in which the rotational DOFs are partially 
'switched on/off'. 

For each of the models considered, the following ex- 
pression: 



V(r) 
kT 



=a 



is fitted within line thickness to the tabulated effective 
potential. 

B. Three-site interaction model 



The coarse-graining procedure described in the paper 
converges after 12 iterations, yielding the effective po- 
tentials shown in Figure [14l 




FIG. 14: Effective pair potentials for the 3-sitc water model. The 
dashed, dashed-dot, and full line curves indicate the O-O, H-H, and 
O-H interactions, respectively. 



C. Two-site interaction model 



The coarse-graining procedure converges after 25 itera- 
tions; the resulting effective potential is shown in Fig- 
ure [T5l 



aie -ln(r-n) +aze -b2(r 



(20) 



D. One-site interaction model 

The coarse-graining procedure converges after 8 itera- 
r2 ^i;ions. The resulting potential for the effective interaction 
between the center-of-mass of two molecules is shown in 
Figure [16] 
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FIG. 15: Effective pair potentials for the 2-site water model. The 
dashed, dashed-dot, and full line curves indicate the O-O, dipolc-dipole, 
and O-dipole interactions, respectively. 
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FIG. 16: The effective interaction between the center-of-mass of two 
molecules as obtained for the single-site water model. 
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